{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "155aacb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "import random\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "39774e7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def trnd(mu, m, n):\n",
    "    ax = np.linspace( -3, 3, 100)\n",
    "    ay = list(stats.t.pdf(ax, mu))\n",
    "    return np.array([random.sample(ay,T), random.sample(ay,T)]).reshape(m, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "a4bc1bb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 10 # Number of observations per simulation\n",
    "K = 10 # Number of resamplings per simulation\n",
    "N = 1000 # Number of simulations\n",
    "# For simplicity, we set var(epsilon) = var(x) = 1\n",
    "betas = np.zeros([4,N])  # Empty matrix to store betas\n",
    "for n in range(N):\n",
    "#     dat = np.random.randn(T,2) # Normal distribution\n",
    "    dat = trnd(4,T,2)  # Student's t-distribution with 4 dof\n",
    "    x = dat[:,0]\n",
    "    eps = dat[:, 1]  # extract x and eps\n",
    "    \n",
    "#     y = x + eps  # true model for y is linear\n",
    "    y = x**2 + eps  # true model for y is quadratic\n",
    "    rhs = np.array([np.ones(T), x]).reshape(T, 2)  # set up rhs matrix for ols\n",
    "    betaols = np.dot(np.dot(rhs.T, rhs)**(-1), np.dot(rhs.T, y))  # estimate via OLS\n",
    "    \n",
    "    # DO THE RESAMPLING PROCEDURE\n",
    "    betas_tmp = np.zeros([2,K])  #  Empty matrix to store betas\n",
    "    for k in range(K):\n",
    "        xy = [[ix, iy] for ix, iy in zip(x, y)]\n",
    "        dat_tmp = random.sample(xy,T)  # resample with replacement\n",
    "        dat_tmp = np.array(dat_tmp)  # shape=(T, 2)\n",
    "        x_tmp = dat_tmp[:, 0]\n",
    "        y_tmp = dat_tmp[:, 1]\n",
    "        rhs = np.array([np.ones(T), x_tmp]).T  # set up rhs matrix for ols\n",
    "        beta_tmp = np.dot(np.dot(rhs.T, rhs)**(-1), np.dot(rhs.T, y_tmp))  # estimate via OLS\n",
    "        betas_tmp[:, k] = beta_tmp  #  write result to matrix\n",
    "\n",
    "    betas_tmp = np.mean(betas_tmp,1)  # average betas across resamplings\n",
    "    betas[:,n] = np.concatenate([betaols, betas_tmp])  # write ols and resampled betas to matrix\n",
    "\n",
    "\n",
    "kurtosis = stats.kurtosis(betas)   # 陡峭度\n",
    "skewness = stats.skew(betas)  # 偏斜度\n",
    "betas_mean = np.mean(betas, 0)\n",
    "betas_var = np.var(betas, 0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "af2e0bfc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.5829041 ,  4.19104428,  2.03822993, 10.25154342])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.concatenate([betaols, betas_tmp])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
